M=15;
point = [7, 8];
Sgrid = zeros(M,M);
for i=1:M
    for j=1:M
        Sgrid(i,j) = exp((-(norm([i,j]-point))^2)/M);
    end
end
surf(Sgrid);
Sgrid = Bivariate_4point(Sgrid,1/16,4);
figure(2);
surf(Sgrid);
